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Abstract 



We derive the first four terms in a series for the order paramater 
(the stationary activity density p) in the supercritical regime of a one- 
dimensional stochastic sandpile; in the two-dimensional case the first 
three terms are reported. We reorganize the pertubation theory for the 
model, recently derived using a path-integral formalism [R. Dickman 
e R. Vidigal, J. Phys. A 35, 7269 (2002)], to obtain an expansion 
for stationary properties. Since the process has a strictly conserved 
particle density p, the Fourier mode N~ l il)k=® - ► P, when N — > oo, 
and so is not a random variable. Isolating this mode, we obtain a 
new effective action leading to an expansion for p in the parameter 
k = 1/(1 + 4p). This requires enumeration and numerical evaluation of 
more than 200 000 diagrams, for which task we develop a computational 
algorithm. Predictions derived from this series are in good accord with 
simulation results. We also discuss the nature of correlation functions 
and one-site reduced densities in the small-K (large-p) limit. 
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I. INTRODUCTION 



Sandpiles are the principal examples of self-organized criticality (SOC) [1-5]. 
Sandpiles with a strictly conserved particle density (so-called fixed-energy sandpiles 
or FES [6]), exhibit an absorbing-state phase transition [7-9], rather than SOC, and 
have recently attracted much interest. Until now, most quantitative results for FES 
have been based on simulations [10-15], an important exception being the solution 
by Priezzhev et al. [16] of a directed, fixed-energy version of the Maslov-Zhang model 
[17], via the Bethe ansatz. Recently, a time-dependent perturbation theory based 
on the path-integral formalism was derived for a stochastic sandpile [18]. In [19] the 
series expansion for the one-dimensional case was extended using operator methods. 

In the present work, the perturbation theory developed in [18] will be reformu- 
lated, leading to an expansion for stationary (t — > oo) properties instead of the short- 
time expansion obtained previously. The expansion parameter is k = 1/(1 + 4p), 
where p denotes the particle density, a conserved quantity. 

Our analysis depends on two basic tools. One is an operator formalism for 
Markov processes, of the kind developed by Doi [20], and which has been applied to 
various models exhibiting nonequilibrium phase transitions [21-25]. The second is 
an exact mapping, devised by Peliti, of a Markov process to a path-integral repre- 
sentation [26,27]. This approach is frequently used to generate the effective action 
corresponding to a process, for subsequent analysis via renormalization group (RG) 
techniques. In the present instance our immediate objective is not a RG analysis 
but an expansion for the order parameter. In the path-integral formalism the prob- 
ability generating function is written in terms of functional integrals over the fields 
ip(x,t) (whose expectation is the particle density at site x), and an auxiliary field 
ijj(x,t). Our reformulation of the effective action is based on the observation that, 
due to particle conservation, the Fourier mode N~ 1 ipk=o is not a random variable, 
but rather has the fixed value p when N, the number of lattice sites, goes to infinity. 

We consider Manna's stochastic sandpile in its fixed-energy (particle-conserving) 
version [12,18,28,29]. The configuration is specified by the occupation number n at 
each site; sites with n > 2 are said to be active, and have a positive rate of toppling. 
When a site topples, it loses exactly two particles ("grains of sand"), which move 
randomly and independently to nearest-neighbor (NN) sites. (Any configuration 
devoid of active sites is absorbing, i.e., no futher evolution of the system is possible 
once such a configuration is reached.) In this work, as in [18,19], we adopt a toppling 
rate of n(n — 1) at a site having n particles, which leads us to define the order 
parameter as p = (n(n — 1)). While this choice of rate represents a slight departure 
from the usual definition (in which all active sites have the same toppling rate), 
it leads to a much simpler evolution operator, and should yield the same scaling 
properties [18]. Preliminary simulation results [30] indicate that in one dimension 
the model exhibits a continuous phase transition at p c = 0.9493. 
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The balance of this article is organized as follows. In Sec. 2, we discuss the 
reorganization of the action, and in Sec. 3 develop a perturbation expansion for 
the activity density in the supercritical regime. Sec. 4 presents the diagrammatic 
expansion rules and the resulting expansion. Predictions for the activity density are 
reported and compared against simulation in Sec. 5, while in Sec. 6 we examine 
correlation functions and higher moments of the density. In Sec. 7 we present a 
brief discussion of our results. 



II. EFFECTIVE ACTION 



As shown in [18], the master equation for the stochastic sandpile can be written 
in the form 



d | g) 
dt 



(1) 



where 



(2) 



{ni} 



where p({rii},t) is the probability of the configuration having occupation numbers 
{rii} and | {rii}) is the direct product of states | rij), representing exactly rij particle 
at site j. In one dimension, the evolution operator takes the form 



-(7Ti_i + 7T i+ i) 2 - 7lf 



(3) 



Here Oj e 7Tj are, respectively, destruction and creation operators associated with site 
i, defined via 



a* | rii) =ni\rii-l) 



(4) 



and 



7Ti | Hi) =| Hi + 1). 



(5) 



As shown in [18], the evolution operator in Fourier representation is given by 



ki,k 2 ,k 3 



(6) 
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with uJk u k 2 — 1 — cos /ci cos/i^; the sums are over the first Brillouin zone. 



m 



As explained in [18], the expectation of any observable A({rii}) can be written 
terms of a functional integral 

(A) = jv^AQ^M, (7) 



where A(ip,ip) is a function of the fields ip and ip corresponding to observable A, 
and 



= exp 



-N' 1 f dt V 4>ki>-k + t did = exp [-iV" 1 f dt C G + f dt d 
Jo ^ Jo 1 Jo Jo 



(8) 



with the interaction given by 



kito,kz 

- 2N~ 2 Wfci,0^fci^fc2^-fci-fca- (9) 
fel,fe 2 



The field ^ is closely related to the occupation number [27]. In particular, the 
activity density is given by 

p{t) = N -ij2( nj ( nj _ 1)) = N- 2 j2M-k) , (10) 

j k 



while for the particle density we have 

<P = N- 1 Y / (n J )=N- 1 (4> k=0 ) . (11) 



In [18] equations (7) - (10) serve as the starting point for a diagrammatic expansion 
of p(t) in powers of time. We now show how these relations may instead be used as 
the basis for an expansion of the stationary activity density = lim t ^ 00 p(t). 

In writing equations (8) and (9) we have assumed a Poisson-product distribution, 
with expectation p, for the initial occupation numbers nj. Thus (ipk=o) — Np, a 
constant of the motion, since the number of particles is conserved. In the infinite- 
size limit, the law of large numbers implies that N~ x il)k=o — P, and is no longer a 
random variable. We may therefore isolate all terms with k = in equation (9) 
setting each factor jV _1 -?/> fc=0 equal to p. (Observe as well that ^pk=o, the variable 
conjugate to ipk=o, is no longer needed.) As a result of this procedure ip] assumes 
the form 
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Q[ijj, ip] = exp — N' 



-i 



Jo 




(12) 



with 



7 fc = 4 p u_ kfi = 4 p (1 - cos fc) 



(13) 



and the modified interaction 




Z UJ k z ,-k 1 -k 2 -k 3 ^k 1 ^k2^k 3 ^-k 1 -k 2 -k 3 



ki,k 2 ,k 3 ^0 



2pN 2 W fc3,-fcl-fc3^fcl^fcs^-fcl-fc3 ~ 2iV 2 Z 



u-k 1 -k 2 ,oipk 1 '*Pk 2 il>-k 1 -k 2 (14) 




P 2jV 1 S U k 3 ,-k 3 ^k 3 ^-k 3 



k 3 ^0 



(15) 



Here it is understood that none of the wavevectors associated with the fields ip and ip 
may be zero. The bilinear part of the action in equation (12) represents independent 
diffusion of particles at rate 4p [27]. The appearance of diffusion at rate 4p in Cq 
may be understood intuitively as follows. The rate of diffusion events at a given 
site is n(n — 1), i.e., twice the number of distinct pairs, so that the diffusion rate 
per pair is 2. The diffusion rate per particle is the twice the diffusion rate per pair 
times the number of pairs per particle, or 4 (n — 1) ~ An ~ Ap if p ^> 1. Unlike 
the original representation of equation (8), the important control parameter p now 
appears explicitly in the action, rather than being "hidden" in the initial probability 
distribution. It is worth noting that this reorganization of the action is not readily 
implemented in the operator representation, equation (3), because in this case it is 
the operator iV -1 J2i ^i^i that assumes a fixed value p. 



III. PERTURBATION EXPANSION 




and the basic contraction or propagator is 



(^ k '( u )Ms))o = N6 k >_ k G(u - s)e 



•7fc(«-s) 



(17) 
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where O represents the step function. As usual in this formalism, 0(0) = [27]. 
The free expectation of n fields if) and n fields ip is given by the sum of all possible 
products of n contractions. 

The expectation of an observable can be written in the form 

(A) = (Aefi c,( *} (18) 

which can be expressed in terms of free expectations if we expand the exponential. 
In this expansion, each field ipkij) must be contracted with a field if)-k(j ), with 
r > t. At n-th order there is a factor of l/n\ and integrations / dti- ■ -dt n over the 
interval [0, £]. We impose the time ordering t > t 1 > t 2 > . . . > t n > 0, thereby- 
cancelling the factor l/n\. We adopt a diagrammatic notation [18] in which fields 
if) (if)) are represented by lines entering (leaving) a vertex. All lines are directed 
to the left, the direction of increasing time. The first term in Cj, equation (14), 
corresponds to a vertex with four lines ( "4- vertex" ) , the second and third to vertices 
with three lines ( "3- vertex" ) , while the fourth, with two lines exiting, will be referred 
to as a "source." figure 1 shows the vertices associated with Cj, as well as the "sink" 
corresponding to the observable p. Vertex b will be called a "bifurcation" and c 
a "junction". In this way, the activity density 

P = N- 2 Y,(M-k) = A^ 2 £(«- fc e/o dt'c' l)o (1Q) 
fc fc 

takes the form 

p = +N- 2 (if)t =0 eJ» dt ' c '') + N- 2 J2(^- k eti dt ' C '')o 
= P 2 + N- 2 £<^_ fc e£ dt ' C >) . (20) 



Consider the first order term. From figure 1 it is evident that the only vertex 
that can be contracted with the sink (without leaving dangling lines) is the source. 
This simple loop, shown as the first diagram on the right hand side of figure 2, makes 
the contribution 

-2p 2 N-^^,-k 

k Jo 

= 1 f ^(i + cosA;)^ 1 -^-!] 

= l{e^[I (8pt) + h(8pt)]-l}. (21) 

where the prefactor 2 is a combinatorial factor and I v denotes the modified Bessel 
function. Here we used 
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Thus this diagram yields the contribution identified in Ref. [18] as p m ax{t), the sum 
of all contributions at order n = 1,2,3, ... proportional to p n+1 , the highest power 
of p allowed at a given order. In the limit t — > oo the contribution to the activity 
from this term is — p/4. 

To study the stationary regime it is convenient to use the Laplace transform. For 
example, the Laplace transform of the contribution due to the simple loop, equation 
(21), is 

2p 2 r n dk 1 - cos 2 k 
s J — 7r 2tt s + 8p(l - cos k) 

where s denotes the transform variable. Using the property lim^oo f(t) = 
lim^o sf(s), we obtain the limiting contribution — p/A directly. 

Consider an arbitrary diagram D of n vertices, and denote the time- dependent 
factors in its contribution to p(t) by /d(0- The Laplace transform of this contribu- 
tion has the form 



f D ( S )= dte~ St dh dt 2 ... dtne -a 1 (t-t 1 )- a2 (t 1 -t 2 )...-a n (t n . 1 -t n ) 

Jo Jo Jo Jo 

roo roo roo 

= (fo ^ / ^j- e -{oci+s)(t-ti)-(a2+s)(ti-t2) — {a„+s)(t n - 1 -t„)-st n 

Jt! Jt 2 JO 

= [s(«i +s)(a 2 + s)... K + s)}-\ (24) 



where the a,- are functions of the wavevectors. Then we have 



f D = hm f D (t) = ft -■ (25) 



The factors a, may be determined via the following procedure. Draw the diagram 
with all vertices in order, and draw vertical lines through each vertex. Then ctj is the 
sum of the factors 7 9 for all propagators between the vertical lines associated with 
vertices % and i — 1 (here t — = to), regardless of whether or not these propagators 
link vertices i and i — 

For example, a diagram composed of n simple loops (see figure 2) makes a con- 
tribution of 

(-l) n 2> 2 
s 

7 



n dk 1 — cos 2 k 



2n s + 8p(l — cos k) 



(26) 



to p(s), and so its contribution to p^ is 

BgV _ (_i)V— (27) 
(8p)« " l i} P (4p)«- 

Summing on n, we find the contribution due to this sequence of diagrams to the 
reduced activity p = lim t _ >00 (p/p 2 ): 

°° /-l\ n I 



In certain cases it is straightforward to replace a simple loop with the infinite sum 
of 1, 2, 3, ... loops. This procedure, illustrated graphically in figure 2, will be called 
dressing a loop. 

Figure 3 shows a three-vertex diagram not included in the sequence equation 
(28). It makes the following contribution to p: 

- 32 P r M (1 + cos k) r ^_ l-cosgcos(k-q) 
Ap(8p) 2 J -it 2tc J -tv 2tt 3 — cos q — cos k — cos(k — q) 



The integral over wavevector q arises frequently in the diagrammatic series and can 
be evaluated in closed form: 



/7T 
-7T 



77 dq 1 — cosgcos(/c — q) 



2n 3 — cos q — cos k — cos(k — q) 



1 

2 



"3-c- v /(l-c)(7-c) 



1 + c 



+ 



/ 1-c 
7- c 



(30) 



where c denotes cos/c. 



In any diagram (beyond the set included in figure 2), we may insert any number 
of loops immediately to the right of the sink. That is, the sink may be replaced by 
a dressed loop. The same applies to the rightmost source, vertex n. The result is 
that the contribution of the original diagram is multiplied by [4p/(l + 4p)] 2 . Once 
this factor is included, no diagram with a 4-vertex immediately to the left of the 
rightmost source (i.e., in position n — 1) or immediately to the right of the sink 
(position 1) need be included in the series. 



IV. DIAGRAMMATIC ANALYSIS 



To begin we define the rules for constructing diagrams in the series for p [18]. 
[Since there is exactly one factor of iV -1 associated with each wavevector sum, all 
of the latter may be changed to integrals, using equation (22).] 
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1. Draw all connected diagrams of n vertices and a sink to the left of all vertices; 
the rightmost vertex must be a source. Each line exiting vertex j must be contracted 
with a line entering some vertex i < j. There is a factor 5k' -k associated with each 
such internal line, where k is the wavevector exiting vertex j and k' the wavevector 
entering vertex i. The requirement that all lines be contracted leads to the condition 
2{n s — 1) + rib — n c = 0, where n s is the number of sources, rib the number of 
bifurcations, and n c the number of junctions. 

2. Each diagram possesses a factor of (—1)" and a combinatorial factor reflecting 
the number of ways of realizing the contractions. In the series for p, this factor is 
given by 2 , with C = 1 + n% + 2n 4 + n s — £, where is the number of 3-vertices 
(of either kind), n 4 the number of 4- vertices, and t the number of simple loops. 

3. Associated with each bifurcation is a factor 2pu kl ^ 2 = 2p[l — cos k\ cos k 2 ] ■ 
Each junction carries a factor 2oUk,o an d each 4- vertex a factor uJk lt k 2 - Each source 
carries a factor of p 2 u)k-k- (The fcj denote the wavevectors exiting the vertex.) 

4. There is a factor f D resulting from the time integrations, as discussed above. 

5. Replace the sink and rightmost source with dressed loops, leading to the 
factor [4p/(l + 4p)] 2 mentioned above, and exclude all diagrams with a 4-vertex in 
position 1 or n — 1. 

6. Integrate over all wavevectors. 

Collecting the factors of p and 1/p associated with the various vertices, f D , and 
the factor of p~ 2 in the definition of p, we find that each diagram in the series for 
p contains an overall factor p~ r where r = n — rib — 2(n s — 1). Using the relation 
2(n s — 1) + rib — n c = 0, we have r = n — n c . 

In order to take advantage of our simple results for the sum of an infinite set of 
diagrams represented by the dressed loops, we adopt k = (l + 4p) _1 as the expansion 
parameter rather than p. Noting that 4p/(l+4p) = 1 — k, and that 1/p = Ak/(1 — k), 
we see that the first order diagram (i.e., the single dressed loop of figure 2) carries a 
factor of 4/(l+4p) = 4/t, while diagrams at higher order carry a factor [4p/ (l+4p)] 2 , 
so that at order l/p r there is an overall factor of (4/t) r /(l — k,) t ~ 2 . Thus diagrams 
oc K r contribute at order r and all higher orders. Diagrams in this class must have 
at least r + 1 vertices and no more than 3r — 2 vertices. 

Enumeration of diagrams at a given order involves (1) identifying all allowable 
sequences of n vertices, and (2) identifying all possible sets of connections between 
vertices, for each sequence. For diagrams with n > 3 (i.e., those not included in the 
simple dressed loop of figure 2), vertex 1 (nearest the sink) must be a junction. (As 
explained above it cannot be a 4-vertex. If it were a bifurcation, the wavevector of 
the single line entering this vertex would of necessity be zero, but such terms have 
been excluded from the action.) For similar reasons, vertex n — 1 must be either a 
source or a bifurcation. Once the vertex sequence has been fixed, all possible sets 
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of contractions of outgoing and incoming lines must be enumerated. The single line 
exiting vertex 1 must, naturally, always terminate at the sink. 

The enumeration of sequences and connections is readily codified in an algorithm 
that may be implemented via computer. In our routine, for each n and r, all 
sequences (subject to the above limitations) are constructed. Then all possible 
connections are generated, by simply running through all termination points for 
each line independently, and rejecting those sets that result in uncontracted lines. 
In this way we are able to enumerate all the diagrams at a given order. 

A diagram is specified in terms of its bond set {(i>i, v[), (v m , v' m )}, where Vj and 
v'j > Vj are the terminal vertices of line j, with j = for the sink. Thus the diagram 
of figure 3 can be written: (01) (12) (12) (23) (03). (The computer algorithm 
was verified against hand enumeration up to third order.) Since the number of 
diagrams grows very rapidly, we extended the routine to perform the wavevector 
integrations for each diagram generated. This entails construction of the numerator 
and denominator of the integrand, which are products of factors involving the cosines 
of various linear combinations of wavevectors. The numerator is a product of factors 
associated with each vertex, as noted in item 3 above. The denominator is a product 
of factors associated with each interval between vertices. These factors are readily 
determined, given the vertex sequence and set of connections. Note that there is 
one free wavevector ki associated with each vertex, except for junctions, so that the 
number of wavevector sums is r. In the latter case the wavevector exiting is equal 
to the sum K of those entering. The lines exiting a source carry wavevectors k! and 
—k', where k' denotes the new associated wavevector. In the case of a bifurcation or 
a 4-vertex we may take the wavevectors of the lines exiting as k' and K — k! where 
K denotes the wavevector entering (or the sum of the wavevectors entering, in the 
case of a 4-vertex). Thus we see that the construction of the integrand (including 
associated numerical factors) is a straightforward task that can also be codified in 
a computational algorithm. The integrals over the ki are evaluated numerically 
using a midpoint method [31]. Based on results for varying number of intervals in 
the numerical integration, we are able to determine the resulting coefficients with a 
relative uncertainty of about 10~ 4 or less. 



V. RESULTS 



We have carried out the expansion for p to order k 4 . Call the number of n- 
vertex diagrams at order k t N n>r and the contribution of this set of diagrams to the 
coefficient of n r /(l — k) t ~ 2 in this series b n y, these values are reported in Table I. 

The diagrammatic expansion yields the following expression for the stationary 
activity density 

Poo = 1 — k — 1.788 040/t 2 - 4.414 481/t 3 - 14.632(2)/t 4 + 0(k 5 ) (31) 



10 



In figure 4 we compare equation (31) and the results of Monte Carlo [19] simulations 
using systems of up to 800 sites. (For each p value, simulations are performed for 
various system sizes and the results extrapolated to the infinite-size limit.) For 
p > 3 the difference between the series expression and simulation is less than 0.1%. 
In Ref. [19], a similar degree of precision is obtained by extrapolating (using Pade 
approximants) a 16-term series (in powers of t) to the infinite-time limit. The present 
series of four terms appears to furnish (without transformation or extrapolation), 
information equivalent to that obtained from a much longer series in powers of t. It 
is worth noting that while the time series is divergent, the present series appears to 
be convergent for small values of n, and it is natural to interpret the first singularity 
on the positive-/? axis as marking the phase transition. 

With a series of only four terms it is of course difficult to draw firm conclusions 
regarding the location of the critical point. We nevertheless analyze the series via 
Pade approximants [32]. The [2,2] approximant is the best behaved and is in excel- 
lent accord with simulation for p > 1.5. It yields a critical value of p c = 0.8677(3). 
(The [1,3] and [3,1] approximants give p c = 0.668 and 0.702, respectively.) It is 
usual to analyse the Pade approximant to the series for the derivative of the log- 
arithm of the order parameter (dlnp/dn in the present instance), as this function 
should exhibit a simple pole at the critical point. The [2,1] approximant does in 
fact give an improved estimate of p c = 0.9069 (about 5% below the value found 
in simulations), while the [1,2] approximant yields p c = 0.860. (The residue at the 
pole of the [2,1] approximant is 0.198, well below any of the numerical estimates for 
the critical exponent f3 that have been reported, which suggest (3 ~ 0.4 [12-14].) 
The series prediction is compared against simulation in figure 4. For p > 3, series 
and simulation agree to within uncertainty, i.e., with a relative error of 10 -4 . We 
note that the four-term series for the stationary activity yields results of accuracy 
comparable to that obtained from the 16-term series in powers of time. The latter, 
when extrapolated to t = oo, gives p c = 0.906 [19]. 

The chief barrier to extending the series is the rapid growth in the cpu time 
required in evaluating the multiple integrals over wavevectors, combined with the 
explosive growth in the number of diagrams. (Enumeration of the diagrams repre- 
sents a small faction of the computing time.) Thus the present approach does not 
appear viable for r greater than four. 

For similar reasons the analysis of the two-dimensional case is restricted to r < 3. 
As explained in Ref. [18], the formalism remains valid in d dimensions if we replace 
all factors uik,k' with 




(32) 



where 



1 



d 



A d (k) 



d 



cosk a . 



(33) 
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Thus 7fc in equation (13) becomes 

7 k = 4p[l-A d (k)] 



(34) 



The expansion involves the same set of diagrams in any dimension; only the inte- 
grals change, with the wave vectors now ranging over the first Brillouin zone in d 
dimensions. 

In two dimensions our result for the stationary activity density is: 

Poo — 1 — k — 1.704 155/t 2 - 3.7292/t 3 + 0(k a ) (35) 

The series prediction is compared against Monte Carlo simulation in figure 5; good 
agreement is observed for p > 1.5. The [2,1] and [1,2] Pade approximants to the 
three-term series for p yield critical values of p c = 0.507 and 0.502 respectively, 
whereas the estimate from simulation is 0.715. 



VI. CORRELATION FUNCTIONS AND PROBABILITY DISTRIBUTION 



Consider the stationary expectation (njUj+e) of the product of occupations at 
sites j and j + I. For i ^ this may be written as [18]: 

C(£) = N- 1 £( W+£ > = N- 2 ]T e M (^. k ) (36) 



and separating the k — term as in equation (20) we find 

C(£) =p 2 + N- 2 £ cos hi {i/> k ij>- k efi dt ' c 'i) (37) 



The second term, an infinite sum of diagrams, defines the connected two-point cor- 
relation function G(|£|). The lowest order contribution comes from the one- vertex 
diagram (simple loop) giving 

G (1 \\£\) = - V - r ^ cosAtf (1 + cosife), (38) 
4 J-tt 2tt 



or G^(l) = —p/8 and G < - 1 - ) (|£|) = for \£\ > 1. When the dressed loop is evaluated 
this becomes G^(l) = —np 2 /8. The correlation for sites separated by greater 
than unit distance is 0(k 2 ) or higher. From this result we can draw the following 
conclusions: (1) The nearest-neighbor correlation is negative for large p; (2) For 
large p correlations decay rapidly in space; (3) As p — > oo, the reduced correlation 
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G(£) = G(i)/p 2 decays to zero (as l/p or faster) so that in this limit the site 
occupancies are independent random variables. 

The stationary expectation of (vpj) m (product of m fields at the same site) is 
related to the m-th factorial moment of the one-site occupation distribution. [For 
m = 2 this is seen explicitly in equation (10).] For m = 3 for example, we have 

(n 3 ) F = N- 1 - l)(n, - 2)) = N~ 3 X> fc ^-fc-fc'> (39) 

j k,k' 

which can be written 

(n 3 ) F = p 3 + 3pN~ 2 Y,(M-k) + N' 3 Yl WkM-k-v) (40) 

The second term equals 3p(p — p 2 ) and so is 0(p 2 ) for large p. The third term must 
be expanded in diagrams in which the sink has three lines entering. The lowest 
order diagram thus involves two vertices, a source and a bifurcation, and is 0(p) 
for large p. We see then that (n 3 ) F = p 3 [l + 0(l/p)\ for large p. The same line 
of reasoning shows that the m-th factorial moment approaches p m as p — > oo. In 
this limit the one-site marginal distribution is therefore Poisson with parameter p, 
and by our previous result on independence, the joint probability distribution is a 
product of such distributions. 

We defer a detailed analysis of correlation functions to future work, and stress 
that the main result of the present section is that in the large-p limit, the probabil- 
ity distribution is a product of identical Poisson distributions at each site, as was 
conjectured in [19]. It is readily seen that this remains valid in d > 2 dimensions. 



VII. DISCUSSION 



We derive a path-integral representation and diagrammatic expansion for the 
stationary activity density in a stochastic sandpile with a conserved particle den- 
sity. Because of conservation, the k = Fourier mode of the particle density (and 
associated field ipk) has a fixed value, rather than being a random variable. This ob- 
servation permits us to reorganize the effective action so that the control parameter 
p appears explicitly, rather than being defined implicitly in the initial condition. The 
bilinear part of the action now describes diffusion at a rate 4p. Because of this, the 
propagator carries an exponential factor, and all time integrations can be realized to 
obtain the limiting (t — > oo) activity directly. The ensuing expansion for involves 
the parameter k = (1 +4p) _1 , i.e., this is a large-p expansion. (As noted in Ref. [19], 
the time series is also most useful for large p values.) We are able to sum certain 
infinite classes of diagrams through the device of "dressed loops." Despite this, the 
number of diagrams to be evaluated at each order grows explosively, so that our final 
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calculational result (the activity series to 0(k a ) is quite modest. The fourth-order 
series agrees very well with simulation in the supercritical regime, and yields (via 
Pade approximation) the critical value p c to within about 10%. A similar favorable 
comparison is seen in the two-dimensional case, although the three-term series fur- 
nishes a poorer estimate for p c . Given these encouraging results, it is reasonable to 
hope that extended series will yield quantitative predictions for critical properties. 
We have also used the reorganized expansion to show that in the large-p limit, the 
sandpile is governed by Poisson-product distribution. Our results strengthen the 
conclusion, until now based on simulation and mean-field-like analyses, that fixed- 
energy sandpiles exhibit a phase transition as the particle density is varied. It is 
of great interest to know if the details of this transition can be analysed using the 
operator and path-integral formalisms. 
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TABLES 



n 


r 




7 


3 


2 


2 


-2.384052 


4 


2 


3 


0.596 013 


4 


3 


4 


5.625 989 


5 


3 


49 


-19.520 916 


6 


3 


180 


11.376 647 


7 


3 


306 


-1.896110 


5 


4 


8 


-13.225188 


6 


4 


311 


135.895 511 


7 


4 


3471 


-311.353 


8 


4 


21961 


256.075 


9 


4 


76 261 


-88.685 


10 


4 


136 404 


11.092 



Table I. Number of diagrams N n . r and coefficient 6 nr in the expansion of the activity. 
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FIGURE CAPTIONS 



Figure 1. Vertices (a - d) in the interaction £ and the sink (e) representing the 
activity density. 

Figure 2. Definition of a "dressed loop" as the sum of one, two, three,... simple 
loops joined at 4-vertices. 

Figure 3. A three-vertex diagram. 

Figure 4. Scaled stationary activity density p versus particle density p in one di- 
mension. Upper curve: series prediction, equation (31); the curve labeled [2,1] is 
obtained by integrating the Pade approximant to the series for dliapdn; points: 
Monte Carlo simulation. Error bars are smaller than the symbols. 

Figure 5. Scaled stationary activity density p versus particle density p in two di- 
mensions. Upper curve: series prediction, equation (35); lower curve: [2,1] Pade 
approximant to the series; points: Monte Carlo simulation. 
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